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1.  Introduction 


Considerable  interest  in  the  gallium  arsenide  permeable  base-transistor 
(PBT)  has  developed  since  it  was  first  announced  by  workers  at  MIT  Lincoln 
Laboratories  in  1979  [1].  Initial  calculations  indicated  that  large  values 
of  the  unity-gain-current  frequency  were  possible  and  could  lead  to  the 
development  of  transistor  amplifiers  and  oscillators  that  would  operate  at 
frequencies  well  into  the  millimeter  wave  frequency  regime.  It  was  also 
suggested  that  LSI/PBT  logic  circuits  with  a  power-delay  product  less  than 
one  fJ  was  possible. 

Since  1979,  considerable  work  has  been  performed  mainly  at  Lincoln 
Laboratories  [2]  and  to  a  lesser  extent  elsewhere  [3].  Workers  have  pre¬ 
dicted  maximum  unity-current-gain  frequencies  as  high  as  80GHz  for  silicon  as 
well  as  operating  frequencies  up  to  30GHz  for  a  silicon  logic  inverter  [3]. 
Additionally,  the  operating  characteristics  have  been  shown  to  be  sensitively 
dependent  on  the  magnitude  of  the  donor  concentrations,  its  shape,  the 
structure  of  the  PBT,  emitter-collecting  spacing,  length  of  the  base  region, 
etc.  [4].  When  one  couples  the  expectation  of  superior  device  performance  with 
its  sensitivity  to  device  parameters,  it  becomes  immediately  clear  that  the 
designer  of  a  high  frequency  PBT  is  faced  with  a  significant  dilemma  in 
choosing  an  optimum  design.  To  compound  these  difficulties,  the  designer  of 
the  PBT  is  also  faced  with  the  dilemma  that  the  principle  design  tool, 
numerical  simulations  of  the  PBT  through  solution  of  the  semiconductor  drift 
and  diffusion  equations  (DDE)  may  not  be  adequate  to  the  task.  The  inade¬ 
quacy  rests  with  the  fact  that  present  PBT  design  trends  are  toward  submicron 
scales  where  the  impact  of  velocity  overshoot  enters.  The  importance  of 
velocity  overshoot  contributions  was  recognized  in  one  of  the  first  PBT 
discussions  where  attempts  to  deal  with  it  were  confined  to  modifications  of 
the  field  dependent  drift  velocity  [2J. 

While  attempts  at  incorporating  velocity  overshoot  effects  into  DDE 
simulations  have  and  will  continue  to  be  useful,  important  design  informa¬ 
tion  is  missing.  The  most  critical  missing  element  lies  at  the  core  of  the 
nonequilibrium  overshoot  phenomena.  Namely,  that  carriers  do  not  attain 
equilibrium  values  of  velocity  instantaneously.  Rather,  there  are  spatial 
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and  temporal  lags  -  features  that  are  absent  from  the  DDE  approach.  To 
overcome  these  deficiencies  newer  and  more  fundamental  approaches  have  been 
developed.  One  approach,  simulation  of  the  PBT  through  solutions  to  moments 
of  the  Boltzmann  transport  equation,  was  successfully  implemented  by  workers 
at  Scientific  Research  Associates,  inc.  (SRA)  during  a  Phase  I  AFOSR/SBIR 
study.  The  purpose  of  this  report  is  to  summarize  the  Phase  I  study. 

The  report  is  divided  into  6  sections.  Section  2  is  a  summary  of  the 
Phase  I  technical  objectives.  Section  3  provides  a  brief  description  of  the 
PBT  and  the  equations  used  in  the  Phase  1  study.  Section  4  contains  a  dis¬ 
cussion  of  the  numerical  procedures.  Section  5  contains  the  results  of  the 
study.  Section  6  contains  the  conclusions  and  recommendations. 

2.  Phase  I  -  Technical  Goals 

There  were  several  goals  of  the  Phase  I  program: 

1.  The  principal  objective  was  to  demonstrate  the  feasibility  of 
adapting  an  existing  numerical  algorithm  for  solving  the  moments  of  the 
Boltzmann  transport  equation  (MBTE)  to  the  study  of  the  gallium  arsenide 
permeable  base  transistor.  The  MBTE  algorithm  incorporates  such  key  high 
speed  and  submicron  features  as  temporal  and  spatial  velocity  overshoot. 

Code  modification  was  successful  and  two  computer  runs  were  performed. 

2.  A  second  objective  was  to  establish  an  effective  program  plan  for 
utilizing  the  MBTE  algorithm  in  the  design  and  physics  analysis  of  the  PBT. 

A  program  plan  was  regarded  as  essential  because  the  costs  associated  with 
MBTE  utilization  are  higher  than  those  associated  with  the  simpler  and  more 
limited  semiconductor  drift  and  diffusion  equations.  The  program  plan 
developed  is  direct  and  involves  two  components. 

i.  The  first  component  involves  choosing  a  specific  design  goal, 

e.g.,  the  design  of  a  60  and  94GHz  PBT  amplifier.  This  component  utilizes 
drift  and  diffusion  equation  algorithms  for  preliminary  selection  of 
material  parameters  (e.g.,  doping  variations),  and  device  parameters 
(e.g.,  base  length,  emitter-collector  spacings)  to  achieve  the  design 
goal. 

ii.  The  second  component  involves  using  the  reduced  set  of  parameter 
variations  as  initial  input  to  the  MBTE  algorithm.  The  MBTE  algorithm 
is  then  utilized  to  simulate  the  operation  of  the  PBT  to  provide 
either  corroboration  or  modifications  of  the  DDE  conclusions. 
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3.  A  third  Phase  I  objective  was  to  enhance  an  ongoing  dialogue  with 
one  of  the  principle  PBT  development  groups ;  that  at  MIT  Lincoln  Labora¬ 
tories.  Workers  at  Lincoln  Laboratories  have  expressed  a  keen  interest  in 
(1)  using  the  results  of  the  MBTE  simulations,  (2)  supporting  part  of  the 
development  necessary  to  bring  an  existing  SRA/DDE  algorithm  on  line  for  PBT 
simulations  and  (3)  utilizing  the  SRA/DDE  algorithm  as  an  adjunct  in  both 
the  interpretation  of  their  experimental  results  and  as  a  guide  to  future 
device  design. 

3.  Permeable  Base  Transistor  and  the  Governing  Moment  Equations 

3a.  PBT  Operation 

Figure  1  displays  a  cutaway  view  of  the  PBT  where  one  notes  the 
presence  of  an  array  of  metallic  fingers  placed  between  two  regions,  labeled 
the  emitter  and  collector  regions.  The  individual  metallic  strips  often  are 
a  composite  tungsten  alloy  and  create  a  Schottky  barrier  at  the  metal/semi¬ 
conductor  interface  region.  The  structure  of  each  region  is  identified  in 
Figure  2,  displaying  an  apparent  similarity  to  that  of  a  three-terminal  field 
effect  transistor.  Indeed,  the  operation  principles  of  the  PBT  and  the  MESFET 
are  similar.  (In  terms  of  nomenclature,  the  following  association  is  made 
between  the  two:  emitter  ■*-*■  source,  collector  drain,  base  •<->■  gate.)  Apart 
from  this  similarity,  there  are  important  design  advantages  of  the  PBT  over 
that  of  the  FET.  First,  substrate  injection  problems  associated  with  the  MESFET 
are  nonexistent  by  design  of  the  PBT  Second,  surface  state  depletion  in  the 
source-gate,  gate-drain  region  of  the  MESFET  is  nonexistent  by  design  of  the 
PBT.  Third,  transport  is  with  the  exception  of  the  base  region,  normal  flow, 
thereby  eliminating  high  current  densities  associated  with  corners  of  the 
planar  FET. 

The  characteristic  electrical  properties  of  the  PBT  are  in  the  zeroth 
order  case  its  dc  current-voltage  characteristics  at  any  combination  of 
emitter-base  collector-base  bias  level.  The  current  flowing  from  the  emitter 
to  the  collector  reflects  the  channel  resistance  between  the  emitter  and 
base,  the  base  and  collector  and  the  enhanced  resistance  within  the  base 
region. 

Under  normally  off  operations  the  PBT  is  designed  such  that  the  deple- 


tion  region  surrounding  the  base  penetrates  to  the  bottom  of  the  channel 
minimizing  the  magnitude  of  current  flow  between  the  emitter  and  collector 
regions.  Current  levels  increase  by  applying  forward  bias  levels  to  the  base 
region.  Figure  3  is  an  illustration  of  the  dc  current-voltage  relation  of  a 
normally  off  PBT,  obtained  from  Reference  5.  As  seen  from  Figure  3,  the  PBT 
begins  to  turn  on  at  a  base  threshold  voltage  of  0.132  volts  with  a  collector 
voltage  of  0.8  volts  [5]. 

3b.  Governing  Equations 

The  operating  characteristics  of  the  PBT  rests  on  solutions  to  the 
governing  transport  equations.  The  governing  equations  utilized  in  this  study 
are  the  first  three  moments  of  the  Boltzmann  transport  equation.  These 
equations  are  derived  for  two  different  species  of  carriers,  in  this  case, 
satellite  valley  and  central  valley  carriers  in  gallium  arsenide.  The  govern¬ 
ing  equations  are  obtained  by  taking  the  collisional  invariant  moments  [6] 
of  the  Boltzmann  transport  equation,  viz.,  the  moments  with  respect  to  the 
mass,  momentum  and  energy  of  the  two  carriers.  This  yields  a  set  of  governing 
equations  which  is  similar  in  form  to  the  equations  utilized  for  two-phase 
fluid  dynamic  flow.  The  governing  equations  reflect  the  conservation  laws 
of  mass,  momentum  and  energy  for  the  two  carriers  and  are  often  referred  to 
as  the  moment  of  the  Boltzmann  transport  equations. 

Using  vector  notation,  the  two  particle  conservation  equations  can  by 
expressed  as 


dn, 

dt 


—  V(n(V,)-n,f1  +n2f2 


(1) 


and 


-  V  (n2V2)  +  n(f(  -  n2f2 


(2) 


where  n^  and  nj 


are  the  satellite  valley  and  central  valley  carrier  number 
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densities  respectively  while  and  are  the  corresponding  velocities. 

and  {2  are  the  corresponding  scattering  integrals  for  mass  conservation 
and  in  general  are  functions  of  the  corresponding  carrier  temperature  and 
momenta  [6].  For  the  purposes  of  this  study,  the  dependence  of  all  scat¬ 
tering  integrals  on  momenta  is  neglected.  Defining  a  total  or  global  number 
density,  n,  by 


and  adding  Eq.  (1)  to  Eq.  (2)  yields  the  equivalent  relationship,  i.e., 
the  global  continuity  equation 

-|^  =  -v[n|v;  +  (n-n|)V2] 

Conservation  of  momentum  for  the  satellite  valley  carrier  can  be  ex¬ 
pressed  as 


-^(n,P,)  *  -V(n,V,^)-  V/J-  V  o-,  -  n,e?-n^f3 

—y 

where  the  momentum,  P^,  and  the  electric  field,  E, are  defined  by 


P.  =  m.Vl 


E  =-V^ 


where  m^  is  the  mass  of  the  satellite  valley  carrier,  e  is  the  electronic 
charge  and  ip  is  the  electric  potential.  The  partial  pressure ,  Pj,  is  re¬ 
lated  to  the  satellite  valley  carrier  temperature,  T^,  and  number  density 
by  the  perfect  gas  relationship 


n,kT, 


where  k  is  Boltzmann’s  constant.  f_  is  the  scattering  integral  for  the 


satellite  valley  carrier  momentum.  The  term  V.o^  represents  the  stress 
forces.  In  this  study,  the  stress  tensor,  a^,  is  approximated  by  the 
Stokes  relationship 


o-,  =  -/1,(VV1  +  7vit-1vv,) 


where  y  is  the  viscosity  associated  with  the  satellite  valley  carriers. 
Substitution  of  Eq.  (6)  and  Eq.  (7)  into  Eq.  (5)  and  dividing 
by  m^  yields  the  final  form  of  the  momentum  equation  for  the  satellite 
valley  carriers,  viz., 

d  — •*  — V'O"  n  — - 

— — { n  V  )  =  -V*(nW  J- — L  — - !+  _J_  eV\//  -  n  Vf  (1( 

di  1  1  1  1  1  m(  m(  m(  1,3 

It  is  observed  that  Eq.  (10)  is  a  vector  equation  and  hence  can  in 
general  be  considered  as  Z  independent  equations  where  SL  is  the  number  of 
relevant  physical  dimensions. 

An  analogous  derivation  for  the  central  valley  carriers  plus  the  identity 
of  Eq.  (3)  yields  the  momentum  equation  for  that  species 

d  r  —1  r  ——1  ^@2  e  — 

^■[(n-n,)V2]=  V  [(n-n,)v2v2]-  —  -  -(n-n,)Vt«4  (1: 


where 


Pz  =  (n-n,)KT2 


cr  =  -u  (VV L  +  VVT-  I-V  VJ 
2  ~2  2  2  ^2 


is  the  mass  of  a  central  valley  carrier  and  Pj  *s  corresponding  viscosity 


coefficient,  f^  is  the  scattering  integral  for  central  valley  carriers  and 
is  assumed  to  be  a  function  of  the  central  valley  carrier  temperature,  only. 

There  are  various  forms  in  which  the  satellite  valley  and  central  valley 
carrier  energy  equations  can  be  described.  This  study  choses  to  cast  the 
energy  equations  in  terms  of  the  satellite  and  central  valley  temperatures, 

T^  and  As  with  the  momentum  equation,  we  start  with  the  satellite 

valley  carrier  energy  equation.  This  equation  can  be  expressed  as  a  balance 
between  the  time  rate  of  change  of  the  total  energy  (internal  plus  kinetic 
energy),  the  convection  of  that  energy,  the  pressure,  stress  and  electrical 
field  work,  the  heat  conduction  and  the  production  and  depletion  of  energy 
due  to  the  change  of  satellite  valley  carriers  to  central  valley  carriers 
and  vice  versa.  Mathematically,  the  above  physical  statement  can  be  ex¬ 
pressed  as 


—  [ni(m(— — 1+  u()j  =  -V- [n(V(  (m|  -L--1  4-  u()j  -  V  (/^V()  —  V •  ( cr  V()  +  n(eV\//  V( 


—  V.  V. 


at 


+  V(?,VT,)-|-k[n1T,(5-(n-n|)V6] 


(14) 


where  is  the  specific  internal  energy  of  the  satellite  valley  carriers, 
k^  is  the  corresponding  thermal  conductivity  and  f,.  and  f^  are  the  scattering 
integrals  for  the  satellite  valley  energy  equation  (which  are  again  assumed 
to  be  functions  of  temperature  only).  To  obtain  the  static  temperature 
version  of  Eq.  (14)  requires  the  elimination  of  the  kinetic  terms 
(V  V  ) 

1  *  1  of  Eq.  (14)  through  the  use  of  the  mechanical  energy  equation. 
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The  mechanical  energy  equation  is  obtained  by  dotting  the  satellite  valley 
carrier  velocity,  V^,  with  the  satellite  valley  carrier  momentum  equation, 
applying  the  vector  identities 


Vl  V("|V|V‘‘7  [("|-!FL)V|] 


vv  — 


(15) 
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-  d  -  d  ,  V.- V.  V.V.  dn. 

Vi'aF(niVi>=  dt  2  +  2  at 


and  the  satellite  valley  particle  conservation  equation,  Eq.  (1) 

This  yields  upon  multiplying  by  m^ 

a  v-  v  r  V-V  _ 

aT(m,n,— >  =  -V  l(m.n|  — )VJ  '  '  Vi  (V  ai)  +  n|eV1VV/ 

+  -^Lm1[n|(2f3-f1)-(n-n1)f2] 

Equation  (17)  is  next  subtracted  from  Eq.  (14)  and  with  the  use  of  the 
vector  identities 


V  ■  (o'j  *  Vj )  =  V,-(7-crj)+  <r(  :  Vv( 


V(/?v,)  =  vrvP,  +  /Jv  v, 


the  internal  energy  version  of  the  energy  equation  is  obtained,  viz., 

^(n|u()  =  V -(n.VjU,)-  /]7-vj  -cr,  :  Vv,+ V-(k,VT,) 

+  m,[n,(2f3-  f,)  -  (n - n,)fj  -  k  [n,T,fs  -  (n -  r.,) T2f6] 


It  is  to  be  noted  that  the  electric  field  term  has  disappeared  in  this  form 
of  the  equation;  in  addition,  the  term  a^rVV^  is  often  called  the  dissipation 
term.  Finally,  to  obtain  the  static  temperature  version  of  the  satellite 
valley  carrier  energy  equation,  the  relationship  between  internal  energy  and 
static  temperature 


V  T  kTi 
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and  the  equation  of  state,  Eq.  (8)  must  be  utilized.  This  yields 

2 

upon  multiplying  by  the  static  temperature  energy  equation. 


£  <".V  -  -V  (n7,T,)-|n,T1V  V, V 

+  3V|-VIm,[n|(2f3-f,)-(n-nJ)f2]  -  n,T,f5  +  (n-n,)T2f4 


(22) 


In  an  analogous  manner  the  energy  equation  can  be  derived  for  the 
central  valley  carrier.  The  results  are 


>\l  ■  -V-[(n-n,>V£Tj-f  (n-n,)TeV-Vt-^o-t:7Vt +|^7-(^T£I 
+  3V2  V2m2[(n-n,)(2f4-f2)  +  n,f,j  +  njT,^-  (n-n,)T2f7 


(23) 


where  f ^  and  fg  are  the  scattering  integrals  for  the  central  vail  energy 
equation  and  is  the  thermal  conductivity. 

The  partial  differential  equations,  Eqs.  (1),  (4),  (10) ,  (11),  (22) 
and  (23)  with  equations  of  state,  Eqs.  (8)  and  (12),  and  the  constitutive 
relations  for  the  stress  tensors,  Eqs.  (9)  and  (13),  constitute  the 
governing  equations  utilized  in  this  study.  Eqs.  (8),  (9),  (12)  and 
(13)  can  be  substituted  into  the  six  partial  differential  equations,  thus 
relating  the  dependent  variables  utilized  in  this  study,  viz.  n^,  n  V^,  V^, 

T^  and  In  addition,  a  governing  equation  must  be  supplied  for  the 

electric  potential  ip  as  this  term  occurs  in  both  the  moment  equations. 

The  electric  potential  can  be  related  to  the  total  number  density  through 
a  Poisson's  equation  of  the  form 


V2^  =  4-<n-No> 


(24) 
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where  Nq  is  the  doping  number  density  (given),  €  is  the  permitivity  and  e 
the  electron  charge.  It  is  to  be  noted  that  the  particle  conservation 
equations,  the  momenta  equations  and  the  energy  equations  all  have  time  (rate) 
terms  while  the  Poisson's  equation  for  electric  potential  has  no  such  term. 

In  addition,  if  we  restrict  ourselves  to  two-dimensional  problems,  it  is 
seen  that  the  governing  partial  differential  equations,  Eq.  (1),  (4), 

(10),  (11),  (22),  (23)  and  (24)  relate  nine  independent  variables  for  nine 
equations  (the  momentum  equations  are  vector  equations  and  in  two  dimension 
are  each  two  separate,  independent  equations).  This  governing  system  of 
equations  is  coupled  and  many  of  the  terms  are  nonlinear.  To  solve  such  a 
system  of  equations  requires  an  initial  distribution  of  each  of  the  inde¬ 
pendent  variables  as  well  as  the  various  coefficients  in  the  governing 
equation,  e.g.  cc^,  etc.  In  addition,  boundary  conditions  are 

needed  to  uniquely  define  the  problem.  The  boundary  conditions  will  depend 
on  the  problem  considered  and  discussion  of  that  matter  will  be  deferred  to 
Section  5  (the  results  section)  where  boundary  conditions  are  specified 
for  a  sample  problem.  An  equation  set  as  complex  as  the  governing  system  of 
equations  (which  will  subsequently  be  referred  to  as  the  moments  of  the 
Boltzmann  transport  equations  (MBTE)  obviously  has  no  closed  form  of  solu¬ 
tion,  except  for  the  most  trivial  of  problems.  Hence,  attempts  to  solve 
this  set  of  equations  with  given  initial  conditions  and  boundary  conditions 
was  accomplished  by  means  of  a  numerical  procedure. 


4.  Solution  of  the  Governing  Equations 


The  numerical  method  used  in  the  solution  of  the  governing  equations, 
which  were  derived  in  the  previous  section,  is  based  on  an  application  of  con¬ 
sistently  split,  linearized,  block  implicit  (LBI)  methods  as  developed  by 
Briley  and  McDonald  [7,8].  LBI  methods  have  been  highly  successful  in  the 
field  of  computational  fluid  dynamics  (CFD)  where  they  have  been  applied  in 
obtaining  solutions  to  a  closely  related  system  of  governing  equations, 
the  Navier-Stokes  equations,  [c.f.  Ref.  7].  Thus,  application  of  such  methods 
to  solution  of  the  moments  of  the  Boltzmann  transport  equation  can  draw  on  a 

vast  amount  of  related  experience  generated  using  LBI  techniques. 

LBI  techniques  center  about  the  use  of  a  formal  linearization  procedure 
in  which  systems  of  coupled  nonlinear  PDE's  in  one  space  dimension  are  reduced 
to  a  system  of  linear  equations,  which  upon  application  of  spatial  differenc¬ 
ing,  may  be  expressed  as  a  block  coupled  matrix  system.  The  resulting  system 
may  then  be  solved  efficiently,  without  iteration,  to  advance  the  solution 
in  time.  Steady  solutions  are  obtained  as  the  long-time  asymptotic  solution. 

The  benefits  of  the  procedure  are  retained  for  multidimensional  problems  through 
application  of  ADI  schemes  in  their  natural  extension  to  block  coupled  systems. 

The  ADI  procedures  reduce  the  multidimensional  system  of  equations,  having 

broad-banded  matrix  structures  to  systems  of  one-dimensional  equations  with 
narrow  block-banded  structures  which  are  solved  efficiently  using  fundamental 
block-elimination  methods. 

Briley  and  McDonald  [7]  considered  the  coupled  system  of  nonlinear, 
time-dependent,  multidimensional  equations  given  by 

a? 

In  Eq.  (25),  <{>  represent  the  vector  of  dependent  variables  ^(n^n^.V^T^tJ/)  , 
H C(f> )  and  S(<j>)  are  nonlinear  functions  of  <f> ,  and  D (<p )  is  a  general,  nonlinear, 
multidimensional,  partial  differential  operator.  Equation  (25)  is  first  time 
differenced  about  tn+8At 


h"+'-h"  .3lDn  +  l  +  Sn+')+<l-g)<Dn+Sn) 


At 


(26) 
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where  At  =  tn+^  -  tn.  The  parameter  3=1  for  a  fully  implicit  scheme  or  3  =  0.5 
for  the  Crank-Nicolson  formulation.  The  implicit  level  nonlinear  operators 
H,D  and  S  are  then  formally  linearized  using  a  Taylor  series  expansion  about 
the  explicit  time  level.  For  example 

h"+i=  -</>v  o (At 2j  (27) 

Eq.  (26)  may  then  be  expressed  at  each  grid  point  in  the  solution  domain 
as  a  matrix  equation  of  the  form 

(A-£AtLn)(<£n  +  ,-<£n)  =  At  (  Dn+  Sn)  (28) 


where 


and 


(29) 


n 


L 


(30) 


As  a  result,  the  nonlinear,  coupled  system  of  PDE's  given  by  Eq.  (25)  has 
been  reduced  to  a  block  coupled,  linear  system  of  temporal  difference  equa¬ 
tions  (Eq.  28)  which,  upon  spatial  differencing,  need  only  be  solved  once 
per  time  step  to  obtain  a  solution.  Additionally,  since  the  linearization 
error  is  at  worst  of  the  same  order  as  the  temporal  discritization  error,  the 
linearization  is  not  expected  to  introduce  significant  inaccuracies. 

Application  of  Eq.  (28)  to  second  order  PDE's  in  one  space  dimension, 
using  standard  three-point  spatial  difference  approximations  requires  the 
solution  of  one  block  tridiagonal  system  per  time  step.  Such  a  system  can 
be  solved  efficiently  using  standard  block  tridiagonal  elimination  procedures. 
However,  application  of  the  LBI  algorithm  given  by  Eqs.  (28-30)  to 
multidimensional  problems  results  in  the  loss  of  this  narrow,  block  banded 
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matrix  structure.  The  discretization  of  the  multidimensional  spatial 
operator  results  in  a  broad-banded  metrix  structure,  which,  if  solved 
by  direct  or  iterative  methods,  can  be  extremely  inefficient.  Such 
observations  led  Briley  and  McDonald  [7,8]  to  develop  consistently 
split  LBI  algorithms  for  multidimensional  problems.  The  splitting 
is  accomplished  by  dividing  the  multidimensional  spatial  operator,  L,  into 
one-dimensional  operators  associated  with  each  coordinate  direction. 


L=L.+L2+L3  OD 

Eg.  (28)  is  then  split  following  the  scalar  ADI  development  of 
Douglas  and  Gunn  [9], 


( A-/3Atl_")  ($*-(}?)  =  At  (Dn+  Sn) 

(32a) 

(A  -y3At  Lg)  =  A  (<£*-<£n) 

(32b) 

(A -£AtLj)  (£***-<£")  =  A (<£**-£") 

(32c) 

Here  A<}>  ,  A<(>  and  A(J>  are  intermediate  solutions  of  Eqs.  (32a-32c)  . 
Again,  if  three-point  operators  are  used  to  approximate  the  spatial 
operators,  L^,  each  of  Eqs.  (32a-32c)  will  be  block-tridiagonal  and  can 
be  efficiently  solved.  The  block  size  and  band  width  are  independent  of 
the  number  of  grid  points,  hence  the  computational  effort  required  to 
solve  the  sequence  varies  linearly  with  the  total  number  of  grid  points 
regardless  of  the  number  of  space  dimensions  considered.  For  two-dimen¬ 
sional  problems,  Eq.  (32c)  is  omitted.  Elimination  of  the  intermediate 
steps  in  Eqs.  (32a-32c)  yields 


(A-/3AIL")  A-,(A-/3AtL2)  A-1  ( A  — /3AtL^)  (<£n  +  '-<£n)  =  At  (Dn  +  Sn) 


(33) 


thus 


,n  +  l  .***-  /f,  ■» 

<f>  =  <£  +  &(A  r)  (34 

The  development  given  above  presents  a  brief  outline  of  the  LBI 
method  used  in  the  present  investigation.  A  more  detailed  development, 
as  well  as  in-depth  discussion  of  LBI  methods,  the  linearization  procedure 
and  related  topics  may  be  found  in  the  article  by  Briley  and  McDonald  [7]. 

Application  of  LBI  procedures  to  Eqs.  (1,4,10,11,22,23  &  24)  in 
two  space  dimensions  give  rise  to  an  A  matrix  of  the  form 


A 


0 

0 

0 

0 

0 

0 

0 

0 

0 


(35) 


This  matrix,  due  to  the  absence  of  a  time  derivative  in  Poisson's  equation, 
is  singular,  thus  Eq.  (28)  cannot  be  split  following  Eq.  (32).  However, 
the  partitioning  of  the  A  matrix  indicated  by  the  dotted  lines  in  Eq.  (35) 
suggests  that  if  the  L  matrix  (the  linearized  D-operator)  could  be  simil¬ 
arly  partitioned,  Poisson's  equation  could  be  decoupled  from  the  rest  of 
the  system.  The  remaining  coupled  equations  could  then  be  solved  by  a 
direct  application  of  the  LBI  method  outlined,  followed  by  the  solution  of 
Poisson's  equation,  completing  a  time  step.  To  accomplish  the  decoupling 
of  Poisson's  equation,  it  is  only  necessary  to  lag  the  electric  field 
terms  (the  potential  gradient)  which  appear  in  the  momentum  equations. 

While  this  formally  reduces  the  accuracy  of  the  temporal  integration  to 


14 


O(At),  it  does  not  adversly  effect  the  stability  of  the  solution  algorithm. 
A  similar  procedure  has  been  employed  by  Kreskovsky  and  Grubin  [10}  in 
solving  the  semiconductor  drift  and  diffusion  equations.  Solution  of 
Poisson's  equation  is  performed  using  a  scalar  ADI  procedure  with  cycled 
acceleration  parameters  [7,  10].  The  overall  solution  algorithm  proceeds 
as  follows: 

1)  Initial  and  boundary  conditions  are  specified  for  all  variables 
throughout  the  computational  domain. 

2)  The  continuity,  momentum  and  energy  equation  are  solved  using 

an  LBI  scheme  to  advance  the  carrier  densities,  velocities  and  temperatures 
from  time  tn  to  time  t°+At. 

3)  Poisson's  equation  is  solved,  using  the  carrier  densities  at 
time  tn+At,  to  obtain  the  advanced  time  potential  distribution. 

4)  Steps  2  and  3  are  repeated  until  a  steady  solution  is  reached  or 
until  the  calculation  is  terminated. 

Finally,  a  few  words  with  regard  to  the  spatial  difference  approxi¬ 
mation  and  computational  mesh  are  in  order.  The  difference  approximations 
used  were  developed  for  arbitrary  grid  spacing  in  either  coordinated 
direction.  Standard,  conservative  three-point  approximations  were  used  to 
assure  that  the  fluxes  of  momentum,  energy  and  current  were  rigorously 
concerned  on  the  finite  difference  mesh.  The  actual  mesh  distribution  can 
be  generated  manually  or  analytic  transformations  such  as  that  of  Oh 
[11],  can  be  used  to  cluster  grid  points  where  higher  resolution  is  required. 
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5.  Results 


The  preceding  analysis  was  incorporated  into  a  computer  program 
aid  the  resultant  program  was  used  to  calculate  the  steady  state  solution 
of  the  above  governing  equations  for  a  gallium  arsenide  permeable  base 
transistor  at  two  separate  base  potentials.  The  results  were  then  com¬ 
pared  with  comparable  results  obtained  by  solving  the  drift  and  diffu¬ 
sion  equations  as  reported  by  Bozler  and  Alley  [2].  A  schematic 
of  the  gallium  arsenide  device  investigated  is  shown  in  Fig.  1.  n  type 

emitters  and  collector  layers  are  used.  The  base  is  patterned  thin  film 
tungsten  while  an  n+  substrate  is  utilized.  Electrons  flow  from  the 
emitter  region  through  the  tungsten  grating  into  the  collector  region. 

The  spacing  between  the  tungsten  bases  is  2000  8  and  the  bases  have  a 
thickness  of  200  8.  The  distance  between  the  emitter  and  collector  is 
1.02pm.  The  tungsten  grating  forms  a  Schottky  barrier  and  with  the 
insulating  gallium  arsenide  controls  the  flow  of  the  electrons  from  the 
emitter  to  the  collector. 

For  the  purpose  of  this  investigation,  an  infinite  bank  of  the 
tungsten  gratings  is  assumed,  hence  planes  of  symmetry  exist  on  the  lines 
passing  through  the  center  of  each  gate  and  at  a  point  half-way  between 
each  gate.  Thus,  it  is  possible  to  perform  the  calculation  in  the  compu¬ 
tational  domain  shown  in  Fig.  4.  The  boundary  conditions  used  in  these 
computations  are  noted  on  Fig.  4.  On  the  planes  of  symmetry  (surfaces 
3)  the  symmetry  conditions  are  imposed,  viz.  that  the  first  derivative 
of  all  variables  except  normal  velocity  are  zero.  The  normal  velocity 
components  are  set  to  zero.  On  surface  1,  the  emitter  contact  surface,  the 
tangential  velocity  components  are  assumed  to  be  zero  as  is  the  central  valley 
normal  velocity  component.  The  second  derivative  of  the  satellite  valley 
normal  velocity  component  is  set  to  zero.  Finally,  the  central  valley 
number  density  and  the  total  number  density  are  set  to  their  initial  values 
(to  be  subsequently  discussed) . 

On  surface  2,  the  collector  contact,  all  variables  are  extrapolated, 
i.e.  second  derivations  are  set  to  zero,  except  for  the  electric  potential 


which  is  specified  at  a  value  V  =1.0  volts.  On  the  surfaces  of  the  em- 

ce 

bedded  gate,  surfaces  A  and  5,  all  normal  components  of  velocity  are  set 
to  zero  (no  penetration  of  the  surface)  while  the  first  derivative  of  the 
tangential  velocity  components  are  set  to  zero  (this  yields  a  slip 
tangential  velocity  on  the  surfaces).  For  the  temperature  the  adiabatic 
condition  is  set  for  both  and  T2,  i.e.  the  first  derivatives  of  temp¬ 
erature  are  set  to  zero.  In  addition,  the  first  derivatives  of  the  number 
densities  n-^  and  n  are  set  to  zero.  The  potential  on  the  embedded  gate 
is  set  equal  to  a  constant  whose  value  is 


(36) 


where  <)>  is  the  Schottky-barrier  height  (set  to  0.8  volts),  E  ,  is  the 
JJN  2 

O  b 

band  gap  energy,  T  is  the  ambient  temperature  (set  to  300  K) ,  Np  is 
the  n-type  impurity  concentration  (set  to  1016/cm3),  n.^  is  the  intrinsic 
electron  concentration  and  V  had  values  of  0  and  0.3  volts  for  case  I 
and  II  respectively. 

Since  this  study  was  limited  to  the  steady  state  solution  of  the 
governing  equations  transient  accuracy  was  not  necessary.  Thus,  any 
initial  condition  that  would  allow  for  the  steady  state  solution  to  be 
obtained  was  acceptable.  The  technique  used  for  this  study  was  to  initially 
set  the  values  of  all  velocity  components  equal  to  zero  (quiescent  flow) 
and  to  set  the  values  of  the  satellite  valley  and  central  valley  tempera¬ 
tures  equal  to  the  reference  temperature,  300°K.  The  global  number  density, 
n,  was  set  equal  to  the  reference  value  and  the  satellite  valley  number 
density  was  calculated  from  the  quiescent  flow  equilibrium  values, 
viz.,  by  using  the  definition  of  total  nondimensional  number  density, 
n,  of  unity  and  neglecting  the  temporal  and  convection  terms  in  the  central 
valley  species  equation  to  obtain 


Initially,  the  electric  potential  was  set  by  solving  Poisson's  equation 

with  zero  space  charge.  For  these  calculations,  the  values  of  the 

viscosity  and  thermal  conductivity  were  assumed  to  be  constant  and 

were  assigned  values  of  6. 1034x10“ 11  ^~cm  and  2. 07x10“ 5 

sec  sec  °k 

respectively.  The  effective  masses  of  the  central  valley  and  satellite 

valley  carriers,  and  m2  were  2.0223xl0~28  gm  and  6. 1034x1. ~29  gm 

respectively.  The  permitivity  €  was  chosen  as  9.8xl0-13  ~0-uj~°-mbs-.  The 

joule-cm 

above  values  were  chosen  as  being  representative  values  for  gallium  arsenide . 

To  demonstrate  the  capability  of  the  computer  code  associated  with  the 
above  analysis,  two  cases  were  run  corresponding  to  base  potentials  of  0.0 
and  0.3  volts,  henceforth  to  be  referred  to  as  Cases  I  and  II  respectively. 

For  Case  I,  two  separate  finite  difference  grid  structures  were  used  to  demon¬ 
strate  grid  independence,  i.e.,  that  the  converged  solutions  for  two  separate 
grid  structures  are  essentially  identical.  For  Case  II  only  one  grid 
structure  was  used.  The  strategy  for  obtaining  converged  solutions  for  both 
cases  was  the  same:  (1)  to  take  a  time  step  as  large  as  possible  (i.e.  as 
large  as  possible  while  maintaining  a  stable  solution)  to  obtain  the  basic 
solution  and  (2)  after  the  basic  solution  is  obtained,  to  take  successively 
smaller  time  steps  to  eliminate  the  remaining  errors  and  thus  to  sharpen 
up  the  solution.  If  one  looks  at  the  Fourier  error  analysis  of  the  governing 
difference  equations,  the  taking  of  the  large  time  step  eliminate  errors 
associated  with  the  low  frequency  components  while  the  taking  of  small  time 
steps  can  be  associated  with  the  elimination  of  errors  associated  with  the 
high-frequency  components.  A  physical  interpretation  is  that  the  taking  of 
the  large  time  step  allows  the  originally  quiescent  solution  to  develop 
its  basic  features.  Although  the  taking  of  such  a  large  time  step  will  not 
give  temporal  accuracy,  this  is,  of  no  concern  as  the  steady  state  converged 
solution  is  the  desired  end.  The  use  of  the  smaller  time  step  allows  the 
finer  features  of  the  solution  to  be  obtained.  The  above  can  be  seen  by 
observing  the  right-hand  sides  of  the  linearized  difference  equations.  (From 
an  iterative  point  of  view,  it  is  the  right-hand  sides  of  the  set  of 
linearized  difference  equations, or  as  it  is  sometime  refererred  to,  as  the 
residual,  that  we  are  attempting  to  drive  to  zero  at  which  time  the  steady 
state  solution  can  be  said  to  be  obtained).  Initially,  the  residuals  are 
zero,  but  the  imposition  of  the  electric  potential  results  in  the  movement  of 
carriers,  thus  resulting  in  the  growth  of  the  residuals.  As  the  basic 


features  of  the  solution  develop,  the  residuals  become  larger  after  which 
the  residuals  will  decrease  to  some  non-negligible  level.  This  corresponds 
to  the  elimination  of  the  low  frequency  error.  If  one  were  to  continue  to 
take  the  same  large  time  steps,  the  residuals  would  cease  decreasing  and  in 
fact  this  criterion  is  used  to  determine  when  one  stops  using  the  large 
time  steps.  At  this  point,  the  time  steps  are  decreased  and  the  residuals 
will  again  start  decreasing.  This  corresponds  to  the  ellimination  of  the 
higher  frequency  errors  and  this  process  is  continued  until  the  residuals 
reach  some  arbitrarily  set  negligible  value.  Usually,  one  can  observe  the  solu¬ 
tion  and  it  will  be  seen  that  as  the  time  steps  are  successively  decreased, 
there  is  a  point  at  which  the  solution  for  all  practical  purposes  ceases  to 
change.  It  is  at  this  point,  that  the  solution  is  deemed  to  be  converged 
and  the  calculation  is  terminated. 

For  Case  I,  two  different  grid  structures  were  utilized.  Case  la 
and  lb  respectively.  Because  of  the  geometry  associated  with  the  permeable 
base  transistor,  a  cartisian  coordinate  system  was  used.  Grid  print  distri¬ 
butions  were  specified  by  the  method  of  Oh  [11].  This  method  allows 
the  grid  points  to  be  specified  by  the  use  of  a  sum  of  error  functions.  The 
advantage  is  that  each  error  function  is  centered  about  some  point  and 
if  the  appropriate  wave  lengths  and  amplitude  constants  are  chosen,  the  effect 
of  that  particular  error  function  will  be  limited  to  a  small  region.  Thus, 
one  has  the  ability  of  concentrating  grid  points  in  desired  regions  while 
using  fewer  grid  points  in  other  regions.  As  used  in  this  study,  relatively 
large  numbers  of  grid  points  were  used  in  regions  where  it  was  expected 
to  have  large  gradients  of  the  dependent  variables.  This  procedure  allows 
for  adequate  resolution  of  these  variables  in  these  regions.  The  first  grid 
distribution  used  for  Case  la  had  21  grid  points  between  the  axis  of  symmetry 
and  51  grid  points  between  the  emitter  and  collector  plates.  The  grid  distri¬ 
bution  utilized  for  Case  lb  used  the  same  21  grid  points  between  the  axis 
of  symmetry  and  63  grid  points  between  the  plates.  For  both  cases  tin  > i id 
points  were  concentrated  in  the  region  of  the  base  (from  the  streamwisi  r 
X-direction  perspective)  and  in  the  region  of  the  top  of  the  base  (fror  the 
transverse  or  Y-direction  perspective).  For  Case  II  the  same  grid  point 
distribution  was  used  as  for  Case  lb. 


Converged  solutions  were  obtained  for  Cases  la  and  lb  in  350  time 
steps.  It  was  noted  that  during  the  convergence  process,  the  satellite 
valley  temperatures  became  as  high  as  8  times  the  reference  temperature 
of  300°K.  However,  for  the  converged  solution,  the  maximum  satellite 
valley  temperature  stabilized  at  approximately  2.5  times  the  reference 
temperature.  It  was  found  that  the  maximum  number  of  time  steps  taken  for 
this  case  were  at  least  an  order  of  magnitude  lower  than  for  the  comparable 
maximum  time  step  allowable  when  the  base  did  not  exist.  Examination 
of  the  reason  for  this  was  seen  to  be  caused  by  the  existence  of  the 
relatively  large  ADI  splitting  error  (see  for  example  ref.  7  for  a  dis¬ 
cussion  of  ADI  splitting  error)  in  the  region  of  the  re-entrant  corner  of 
the  base.  Since  similar  phenomena  have  been  observed  in  fluid  dynamics 
calculations  with  re-entrant  corners,  this  phenomenon  was  not  unexpected  [7  ]. 
Since  the  splitting  error  is  proportional  to  the  time  step  squared. 

At2,  the  splitting  error  can  usually  be  controlled  by  taking  a  smaller 
time  step.  This  was  the  technique  used  in  the  Phase  I  study.  The  penalty 
to  be  paid  in  taking  small  time  steps  is  the  use  of  more  time  steps  (and 
hence  more  computer  time)  to  obtain  a  converged  solution.  For  case  lb, 
the  CPU  time  (on  the  University  of  Minnesota  Cray-1  computer)  is  5.4xl0~3 
seconds  per  grid  point  per  time  step  for  nine  equations.  For  case  lb  this 
corresponds  to  approximately  2,500  seconds  (approximately  42  minutes) 
of  CPU  time. 

For  the  Case  II  computations,  the  converged  solution  for  Case  lb  was 
used  as  an  initial  condition.  As  in  case  lb,  a  converged  solution  was  ob¬ 
tained  for  Case  II  in  350  time  steps.  Again,  the  splitting  error  necessi¬ 
tated  the  taking  of  a  smaller  time  step  than  would  have  been  necessary  for 
the  case  without  the  base  region.  Since  the  gradients  are  larger  in  the  re¬ 
gion  of  the  base  when  the  base  voltage  is  increased,  this  case  had  corres¬ 
pondingly  larger  splitting  error  for  the  same  time  step  as  was  used  for  Case 
lb.  Thus,  an  even  smaller  time  step  had  to  be  used. 

The  results  of  the  calculation  are  displayed  as  contour  plots.  For 
the  case  of  0.0  volts  on  the  base  contact,  the  potential  contours  are  shown  in 
Figure  5.  It  is  noteworthy  that  the  largest  potential  gradient  along  the 
line  of  symmetry  within  the  channel  occurs  immediately  downstream  from  the 
base  contact.  The  highest  electric  field  within  this  region  is  approximately 
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40kV/cm.  Within  the  framework  of  the  drift  and  diffusion  equations  this 
would  tend  to  fix  the  mean  carrier  velocity  at  1.0xl07cm/sec  and  assure 
that  well  over  90%  of  the  T  valley  carriers  had  transferred  to  the  satellite 
valleys.  The  percentage  of  T  valley  carriers  along  the  line  of  symmetry 
arising  from  solutions  to  the  moment  equations  is  displayed  in  Figure  6 . 

It  is  apparent  that  considerably  less  transfer  has  occurred  than  that  asso¬ 
ciated  with  the  drift  and  diffusion  approximation.  Indeed,  the  mean  carrier 
velocities  within  the  channel  are  higher  than  those  associated  with  the 
equilibrium  field  dependent  drift  velocity,  and  at  certain  points  exceed 
4xl07cm/sec. ,  Fig.  7.  The  density  contours  for  this  case  are  shown  in 
Fig.  8  and  these  contours,  as  well  as  those  of  Fig.  5  are  noteworthy 
for  their  similarity  to  those  of  Reference  2.  It  may  be  stated  that 
while  the  carrier  and  potential  contours  are  similar,  the  distributions 
of  carriers  between  the  central  and  satellite  valleys  are  different 
from  those  using  the  drift  and  diffusion  equations. 

While  it  is  unlikely  that  a  full  range  of  dc  current  -  voltage  charac¬ 
teristics  will  yield  similar  results  using  the  moment  equation  algorithms , 
and  the  drift  and  diffusion  algorithms,  if  they  do,  the  question  arises  as 
to  the  usefullness  of  extracting  high  frequency  device  performance  charac¬ 
teristics  from  dc  characteristics. 

The  results  for  the  forward  bias  calculations  are  shown  in  Figs. 

9-11.  For  this  case,  the  field  distribution  downstream  from  the  base 
contact  is  somewhat  lower  than  that  of  Figs.  5  and  8  and  the  magnitude 
of  electron  transfer  is  muted  compared  to  zero  applied  base  contacts.  The 
distribution  of  carrier  density  along  the  channel  is  displayed  in  Fig.  10 
and  the  remarks  concerning  the  comparative  results  obtained  with  the  drift 
and  diffusion  code  are  the  same  as  that  associated  with  Fig.  6. 

6.  Conclusions  and  Recommendations 

The  analysis  performed  under  the  Phase  I,  AFOST/SBIR  study  demonstrates 
that  a  recently  developed  algorithm  can  be  successfully  aDDlied  to  ex¬ 
amining  the  electrical  characteristics  of  the  gallium  arsenide  permeable 
base  transistor.  The  study  indicates  that  there  is  a  possibility  for  the 
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dc  characteristics  as  obtained  by  drift  and  diffusion  equation  algorithms 
to  be  similar  to  those  obtained  using  the  moment  equation  algorithms. 

However,  the  carrier  and  velocity  distribution  will  be  different,  and  pre¬ 
dictions  based  on  the  drift  and  diffusion  algorithms  are  likely  to  be  of 
limited  value  for  high  frequency  operation. 

While  the  Phase  I  study  demonstrated  the  feasibility  of  using  the  moment 
equation  algorithm  for  PBT  study,  additional  algorithm  development  is  re¬ 
quired  to  improve  convergence  rates,  particularly  in  the  vicinity  of  the  base. 
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FIGURE  4.  COMPUTATIONAL  OOMAIN  WITH  BOUNDARY  CONDITIONS 
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FIGURE  12.  CONTOURS  OF  CARRIER  DENSITY  -  BASE  0.3  volts 


